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Mammalian skull heterochrony reveals modular 
evolution and a link between cranial development 
and brain size 
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The multiple skeletal components of the skull originate asynchronously and their develop- 
mental schedule varies across amniotes. Here we present the embryonic ossification 
sequence of 134 species, covering all major groups of mammals and their close relatives. This 
comprehensive data set allows reconstruction of the heterochronic and modular evolution of 
the skull and the condition of the last common ancestor of mammals. We show that the mode 
of ossification (dermal or endochondral) unites bones into integrated evolutionary modules of 
heterochronic changes and imposes evolutionary constraints on cranial heterochrony. How- 
ever, some skull-roof bones, such as the supraoccipital, exhibit evolutionary degrees of 
freedom in these constraints. Ossification timing of the neurocranium was considerably 
accelerated during the origin of mammals. Furthermore, association between developmental 
timing of the supraoccipital and brain size was identified among amniotes. We argue that 
cranial heterochrony in mammals has occurred in concert with encephalization but within 
a conserved modular organization. 
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The mammalian cranium displays great morphological 
disparity, reflecting the wide spectrum of ecological 
diversification in the groups A major aspect of skull 
diversification is heterochrony, or changes in developmental 
timing, for which genetic mechanisms are now increasingly being 
understood^'^. Simple alterations in the onset, duration and tempo 
of development are regarded as causes of profound morphological 
changes'*. Until recently, most mammalian heterochronic studies 
have focused on postnatal life, and our knowledge of fetal 
development has largely been restricted to model organisms^. The 
critical stages for examination of organogenesis are fetal or around 
the time of birth^'^, and thus non-model organisms are rarely 
available and difficult to sample^'^. Recent progress in high- 
resolution imaging techniques has provided new avenues to non- 
destructive investigation of fetal and neonatal specimens from 
museum collections. Microtomographic imaging allows 
documentation of the onset of individual bone ossification, a 
powerful marker for tracing perinatal anatomy. This approach has 
served to identify fijndamental differences in postcranial 
osteogenesis between marsupials and placentals*^, and 
unsuspected variation in placental development* ^ However, the 
ancestral patterns of craniogenesis timing and factors behind the 
cranial heterochrony remain largely unknown. In this study, we 
show that the timing of bone formation in the mammalian skull is 
greatly influenced by two factors such as brain size and 
developmental modularity caused by the mode of ossification. 

Encephalization is a central phenomenon in mammalian 
evolution, one that has led to the largest brained vertebrates*^, 
as best exemplified by primates and cetaceans. During the 
early evolution of mammals in the Jurassic, brain expansion 
was associated with the acquisition of the neocortex and 
diversification of sensory faculties*^. Ontogenetically and 
evolutionarily, the expansion of the cranial vault reflects brain 
size increase* ^~*^, as shown also for humans*^. Given such 
somatic integration between the skull and the brain, we tested 
whether the heterochronic changes in embryonic ossification 
reflect the evolution of brain size. Modularity, referring to the 
strong internal integration and weak interactions among 
morphological subsets*^, is another aspect of patterns of 
heterochrony to consider*^. It has been suggested that genetic 
modularity affects the evolutionary dynamics of species, which in 
turn influence the evolution of molecular networks regulating 
morphogenesis*^'*^. However, the link between heterochrony and 
modularity in macroevolution remains largely unknown^^" . 

A comprehensive sampling of museum collections across the 
world using non- destructive micro -computed tomography tech- 
nique produced skeletal developmental sequences for 21 cranial 
elements of 102 mammalian species and 32 non-mammalian 
amniote species (sauropsids). Covering almost all major mam- 
malian groups, our exceptionally large data set was used to 
reconstruct the developmental sequence of the common ancestor 
of mammals, and provide insights into evolutionary patterns of 
skeletal development. Herein, we demonstrate that cranial 
heterochrony reflects the encephalization history of mammals 
and conserved modular organization of skull elements. 

Results 

Reconstruction of ancestral ossification sequence. The ancestral 
conditions and heterochronic changes of ossification sequence 
at all nodes were reconstructed from developmental sequences 
of 134 amniote species (Supplementary Data 1 and 2; and 
Supplementary Table 1) using squared- change parsimony under a 
Brownian motion model of character evolution . The ossification 
sequence inferred for the ancestor of Mammalia is as follows: 
(1) premaxilla and maxilla, (2) dentary, (3) palatine, (4) frontal 



and squamosal, (5) pterygoid, (6) jugal, (7) parietal, (8) nasal 
and ecto tympanic, (9) vomer, (10) exoccipital and goniale, 
(11) basioccipital, (12) lacrimal, supraoccipital and alisphenoid, 
(13) basisphenoid, (14) orbitosphenoid and (15) petrosal (Fig. 1 
and Supplementary Table 2). Inferred heterochronies for all other 
nodes are given in Supplementary Figs 1-21. We also conducted 
an alternative reconstruction by Parsimov-based genetic inference 
(PGi)^"*. This approach treats the sequence as one single, 
complex character and uses the Parsimov algorithm^^ as an 
edit-cost function to optimize ancestral states and sequence 
heterochronies. The inferred sequence for the ancestor of 
Mammalia is: (1) premaxilla, (2) maxilla, dentary, nasal, jugal, 
frontal, parietal, squamosal, vomer, palatine and ectotympanic, 
(3), lacrimal, (4) basioccipital and supraoccipital, (5) pterygoid, 
(6) basisphenoid, (7) orbitosphenoid, alisphenoid and exoccipital, 
8) goniale and (9) petrosal. The PGi sequence was less resolved 
(that is, involving more tied ranks) than that generated 
by squared- change parsimony, but both were mostly similar 
(Spearman's rank correlation analysis, rs = 0.80, n = 21, 
P< 0.0001). Inferred heterochronies for higher taxonomic levels 
are given in Fig. 2, and those for more inclusive nodes are given in 
Supplementary Figs 22-29. The results obtained by squared- 
change parsimony showed that the last common ancestor 
of Mammalia had a more accelerated onset of ossification 
of the vomer, frontal, parietal, basioccipital, exoccipital, 
and supraoccipital compared to non-mammalian amniotes 
(Supplementary Table 2). The PGi analysis showed that the last 
common ancestor of Mammalia had a more accelerated onset of 
ossification of the frontal, parietal, basisphenoid, basioccipital and 
supraoccipital (Fig. 2). 

Ossification patterns and encephalization. Results of correlation 
analysis of the relative timing of cranial ossification (scaled from 0 
to 1) and encephalization quotient (EQ), which is the residual of 
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Figure 1 | Reconstructed ossification sequence of the hypothetical 
common ancestor of Mammalia using squared-change parsimony 
under a Brownian motion model. The skull of the mammaliaform 
Morganucodon^'^'^'^ is used to show the adult bone topology. The skull is in 
lateral view and the lower jaw in lingual view. Septomaxilla, coronoid and 
articular (nnalleus) were not applicable or studied herein. Abbreviations: 
as, alisphenoid (epipterygoid); bo, basioccipital; bs, basisphenoid; 
de, dentary; eo, exoccipital; et, ectotympanic (angular); fr, frontal- 
go, goniale (prearticular); ju, jugal; la, lacrinnal; mx, maxilla; na, nasal; 
OS, orbitosphenoid; pa, parietal; pal, palatine; pe, petrosal; pg, pterygoid; 
pm, premaxilla, so, supraoccipital; sq, squamosal. This material is 
reproduced and modified with permission of John Wiley & Sons, Inc. 
Copyright 1981 (John Wiley & Sons, Inc). 
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Figure 2 | Heterochronic shifts in the onset of skull bone ossification recovered by the Parsimov-based genetic inference (PGi) analysis in 
amniotes. Significant shifts detected in derived nodes connpared with ancestral nodes are sunnnnarized. A, acceleration; D, delay. Numbers in the tree 
represent the detected cranial elennents. 
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an allometric regression of brain weight against body weig 
are given in Table 1. We found that the ossification onset of the 
supraoccipital bone occurs earUer in taxa with higher EQ in 
mammals (Pearson's product moment correlation analysis, 
r= -0.65, n = 48, P<0.0001) (Fig. 3). Similarly, among non- 
mammalian amniotes, a tight correlation was found between the 
developmental timing of the supraoccipital and EQ (Pearson's 
product moment correlation analysis, r = — 0.95, n = 9, 
P< 0.001). These correlations were similarly significant in the 
phylogenetically controlled correlation analysis (Fig. 4; Table 1). 
The results by squared- change parsimony revealed that the 
timing of the supraoccipital development was accelerated in 
ancestral mammals when splitting from non-mammalian 
amniotes (Fig. 5 and Supplementary Table 2). The 
developmental timing of the supraoccipital was shown to be 



significantly earlier in mammals than in non- mammalian 
amniotes ((7- test on phylogenetic independent contrasts^^, 
mammals n — 79, non-mammalian amniotes n = 30, P< 0.001). 



Modularity analysis. Depending on the anatomical identity of 
skull bone elements, the cranial region can be divided into dif- 
ferent modules. Developmentally, the skull can be divided into 
mesoderm and neural crest cell-derived elements^^. Similarly, 
skull bones are classified into either dermal bones or 
endochondral bones, depending on their mode of ossification^^. 
Morphometric analyses of adult mammalian skulls have 
previously identified five phenotypic variational modules: oral, 
zygomatic, nasal, cranial base and cranial vault^^. We considered 
these divisions as hypothetical modules, and tested whether these 



Table 1 | Pearson's correlation coefficients and P-values for 
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Both raw comparisons and phylogenetically corrected comparisons (Felsenstein's independent contrasts) are given. Significance level was set as P< 0.05/21 after Bonferroni correction, and statistically 
significant values are given in bold. Ectotympanic, goniale and alisphenoid of mammals were homologized to angular, prearticular and epipterygoid of non-mammalian amniotes, respectively. Values for 
orbitosphenoid of non-mammalian amniotes are not available because of uncertain homologies of this bone. Values for alisphenoid (epipterygoid) of non-mammalian amniotes are not available, as both 
the ossification timing of this bone and EQ is reported only for one species (Locerto ogilis). 
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Figure 3 | Relation between encephalization quotient (EQ) and supraoccipital ossification timing. Supraoccipital tinning was obtained by 
calculating the relative score of the supraoccipital within all cranial bones. Significant negative correlation was found between supraoccipital timing 
andEQ(r= - 0.65, n = 48, P< 0.0001). 
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Figure 4 | Phylogenetic independent contrasts between EQ and supraoccipital timing. Supraoccipital timing is negatively correlated with the 
EQ (r= -0.60, n = 47, P<0.0001). 



modules are identifiable in the patterns of skeletal heterochrony. 
Our results on pooled species demonstrated that timing of 
ossification of dermal bones is skewed towards earlier 
developmental stages than that of endochondral bones ((7- test, 
n = 1,470, P<0.05) (Fig. 6a). In addition, neighbour- joining 
cluster analysis showed that the skull bones form two evident 
clusters, one cluster consisting explicitly of dermal bones and the 
other of endochondral bones (n = 1,470) (Fig. 6b,c). 



Discussion 

The common ancestor of Mammalia was found to have an 
accelerated onset of ossification of the cranial bones associated 
with the braincase (frontal, parietal, basioccipital and supra- 
occipital) when compared with non-mammalian amniotes. 
Morganucodon (Fig. 1), one of the basal-most mammaliaforms 
(the clade that includes mammals and their closest relatives), had 
a greatly expanded olfactory bulb, olfactory cortex, neocortex and 
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Figure 5 | Reconstructed heterochronic shifts of supraoccipitai deveiopmentai timing. 

under a Brownian motion nnodel. 



Values were reconstructed using squared-change parsinnony 



cerebellum in comparison with non-mammaliaform cynodonts . 
Such a morphological transformation is recognized as the first 
evolutionary pulse of brain expansion in mammalian evolution. 
The expanded brain regions are covered by accelerated bones of 
the skull- roof, suggesting that increased encephalization led 
to quantifiable developmental changes in the skull. Among 
mammals and non-mammalian amniotes, the ossification onset 

6 



of the supraoccipitai bone, which covers the occipital lobe of the 
cerebrum and cerebellum, is more accelerated in taxa with higher 
EQ (Figs 3 and 4; and Table 1). This indicates that the 
developmental timing of the supraoccipitai can predict brain 
size. The developmental timing of the supraoccipitai is more 
precocious, on an average, in mammals than in non-mammalian 
amniotes. Furthermore, our squared-change parsimony analysis 
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Figure 6 | Variation of deveiopmentai timing of ail bones. Results are shown for highly resolved species with more than four ranks, (a) Boxplot 
comparison of the range of variation of relative timing of cranial bones (dermal bones in red and endochondral bones in blue). The 25-75th percentiles are 
shown using a box. The whiskers show the top of the box up to the largest data point <1.5 times the box height from the box. Values further than 1.5 times 
the box height from the box are shown as circles, and those further than three times are shown as stars. Relative ossification timing is indicated to be 
significantly different between pooled dermal bones and pooled endochondral bones (Mann-Whitney U-test, P< 0.00001). (b) Neighbour-joining cluster 
analysis of developmental timing. Bootstrap values were obtained through 10,000 permutations. Cranial bones form two clusters such as a dermal bone 
cluster and an endochondral bone cluster, (c) Topology of dermal and endochondral bones (prenatal Bos taurus). Abbreviations: as, alisphenoid; bo, 
basioccipital; bs, basisphenoid; de, dentary; eo, exoccipital; et, ectotympanic; fr, frontal; go, goniale; ju, jugal; la, lacrimal; mx, maxilla; na, nasal; os, 
orbitosphenoid; pa, parietal; pg, pterygoid; pm, premaxilla, so, supraoccipital; sq, squamosal. 



detected that the timing of the supraoccipital development was 
considerably accelerated in ancestral mammals when splitting 
from non-mammalian amniotes, but then remained constant in 
monotremes, marsupials and the last common ancestor of 
placentals (Fig. 5). Then, its timing was further accelerated in 
multiple lineages independently, including primatomorphans, 
cetaceans, talpids and dipodid rodents (Fig. 5), all of which are 
known as encephalized^^' . Within primates, humans exhibit the 
highest EQ and the most accelerated case of supraoccipital 
development. Among rodents, the supraoccipital is most 
accelerated in jerboas (Jaculus) that possesses the highest EQ 
among the studied rodents (Supplementary Table 3). On the 
other hand, the ossification timing of another major skull- roof 
bone, the parietal, is not significantly correlated with EQ. 
However, it is worth noting that the onset of parietal 
ossification was also accelerated at the common ancestor of 
placentals (Supplementary Fig. 8). Moreover, the parietal is one of 
the earliest bones to develop even in other non-placental 
amniotes, and its ossification timing varies little among 
placentals. Parietal timing possibly reached a plateau at the 
placental ancestor, and therefore it is not correlated with EQ. 

Recent genetic studies have shown that the development of the 
supraoccipital and brain are genetically integrated. The apparent 
link between supraoccipital development and brain expansion 
may be because of the pleiotropic effect of Lmxlb and Dlx5. 



The supraoccipital and interparietal are either absent or severely 
reduced in Lmxlb knockout mice^^. Furthermore, this gene is 
critically required for mid/hindbrain development^^. Dlx5 is 
essential for axonogenesis and nervous system development and 
is reported to be related to Down Syndrome in humans^^. It also 
affects the timing of supraoccipital ossification^"*, and more 
importantly, Dlx5 null mutants explicitly lack the supraoccipital 
and interparietal^^. 

Our modularity analysis demonstrates that timing of ossifica- 
tion of dermal bones is constrained towards earlier developmental 
stages, whereas that of endochondral bones occurs later (Fig. 6a). 
Furthermore, the skull bones form two separate modules, one 
consisting explicitly of dermal bones and the other of endo- 
chondral bones (Fig. 6b,c). It is possible that when sequence 
heterochrony occurs during evolution, developmental timing of 
bones of identical developmental modes are likely to covary, and 
that bones of different developmental modes tend to be more 
independent from each other. On the other hand, neither 
mesoderm versus neural crest origin nor phenotypic modularity 
identified based on adult metric traits^^ appears to be related to 
cranial ossification heterochrony. In the early stages of vertebrate 
development. Hedgehog signalling critically controls the 
differentiation of osteoblasts and the onset of osteoclast activity 
in endochondral bones, while alteration of this signalling has little 
effect on dermal bone formation^^. Such a finding implies the 
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genetic independence of endochondral bones and dermal bones, 
that is, timing of osteogenesis of endochondral bones and dermal 
bones are controlled by somewhat independent gene regulatory 
networks. Together, we suggest that such genetic integration 
constrains cranial ossification timing both ontogenetically and 
evolutionarily. 

Our study highlights the conserved modular organization 
imposed on cranial heterochrony and the evolutionary degrees of 
freedom in this integrated system. Ossification modes fundamen- 
tally constrain the evolvability of cranial development. Although 
the integration of all other endochondral bones is evident, the 
supraoccipital appears to be rather independent from the rest 
(Fig. 6b). This bone exhibits the most variable ossification timing 
among endochondral bones (Supplementary Table 4) and does 
not form a tight ossification timing cluster with other occipital 
elements (that is, exoccipital and basioccipital) (Fig. 6a,b), despite 
the shared somite derivation of all three occipital components'^. 
We suggest that this relative independence of the supraoccipital 
may be because of its tight link with the brain. 

Methods 

Specimen collections. Specimens sampled are held at the Anthropological 
Institute and Museum of University of Ziirich (AIMUZ), Botanical Gardens 
Museum of Hokkaido University (BGHU), Institute of Ecology and Biological 
Resource of Vietnamese Academy of Science and Technology (lEBR), Japan 
Monkey Center (JMC), Kyoto University Museum (KUM), Natural History 
Museum Bern (NMB), Natural History Museum Wien (NMW), Swedish Museum 
of Natural History Stockholm (NRS), National Museum of Nature and Science 
Tokyo (NSMT), Palaeontological Institute and Museum of University of Ziirich 
(PIMUZ), Wildlife Laboratory at Tokyo University of Agriculture (TUA), 
University Museum of University of Tokyo (UMUT), and Berlin Museum of 
Natural History (ZMB). Specimens used in this study are summarized in 
Supplementary Data 1. 

Data acquisition. Ossification sequence data of 21 cranial elements were 
documented. Ectotympanic, goniale and alisphenoid of mammals were homo- 
logized to angular, prearticular and epipterygoid of non-mammalian amniotes, 
respectively-^^. The appearance of bones was assessed non-invasively by acquiring 
shadow images taken by |iCT at the University Museum, University of Tokyo 
(TXS225-ACTIS, TESCO, Tokyo) and at the Anthropological Institute, University 
of Ziirich (|aCT80, Scanco Medical, Bassersdorf). Three-dimensional visualization 
and analysis of shadow images were conducted in Amira 5.3 (Visage Imaging 
GmbH, Berlin, Germany). Supplementary Data 2 lists the acquired sequences and 
those obtained from the literature. 

Phylogenetic framework. The topology was arranged in Mesquite^^. Phylogenetic 
framework of species studied and divergence time are based on molecular evidence 
(Supplementary Table 1). Divergence time has been estimated for all major 
sauropsid and mammalian clades. Therefore, although not completely consistent 
internally, the TimeTree of Life-project"^^ resamples the most comprehensive 
synopsis of molecular-based phylogenetic studies to date. Several authors published 
divergence times of lower taxonomic levels. Those were brought into relation with 
the TimeTree of Life. Therefore, the deepest, most overlapping phylogenetic node 
between the TimeTree of Life-project and the specific study were compared among 
each other and brought into relation. The resulting factor was then used to 
normalize the divergence times of lower taxonomic levels in the specific study. 
Usually these were known only for one or two subclades. Such normalization was 
also performed for data of those chapters in the TimeTree of Life'^^, which show 
inconsistency towards the higher phylogenetic levels of other chapters. Only few 
studies exist that present molecular-based divergence times of the mammalian or 
sauropsid subgroups. Moreover, those studies often do not show nodes that overlap 
with nodes in the phylogeny of the TimeTree of Life or the subclades of our 
taxonomic sampling are not represented. In those cases, the branch lengths 
between the nodes (of unknown age) within a major clade (of known divergence 
time) were evenly distributed (Supplementary Table 1). Based on this strategy of 
dating the composite phylogeny, the significance of our results is particularly high 
on the higher taxonomic levels. Only those are discussed in the present 
contribution. 

Heterochrony analysis. We used two methods, squared- change parsimony^^'^^ 
and PGi^'^, to reconstruct heterochronic changes in amniotes. In the former 
approach, the sequence of each bone is divided by the maximum rank, resulting 
in intervals that are standardized between 0 and 1. Then, squared-change 
parsimony^^ based on a Brownian motion model of character evolution and 



Felsenstein's independent contrasts^^ is used to reconstruct the heterochronic 
changes at all nodes. The analysis was conducted with the PDAP module of 
Mesquite-^^. Divergence times derived from molecular dating were used as branch 
lengths. As the resolution of the sequence can bias the results in this approach, 
well-resolved species with more than three ranks were included. The alternative 
PGi examines the sequence as one single, complex character and uses the Parsimov 
algorithm as an edit-cost function to optimize ancestral states and sequence 
heterochronies. The PGi algorithm computes the lowest cost assignment of the 
ancestral sequences in a two-step, dynamic programming procedure^^. The 
advantage of this approach is that no assumptions are made of the data, outside of 
those made when evaluating the hypothetical solutions^'^. The parameters used for 
the analysis were as follows: 100 cycles, 100 replicates and 100 sequences retained 
at each node. Semi- exhaustive search with 10,000 permutations was performed. 
Such runs were conducted four times independently, and the shortest tree was 
treated as the conservative reconstruction. As the phylogenetic position of turtles is 
still disputed, and as results by PGi can be affected by polytomies, turtles were 
excluded from this analysis. The analysis was conducted using 'ape', 'el 071', and 
'PGi' packages in R^^. 

Comparisons with brain size. We compared the relative timing of cranial 
ossification (scaled from 0 to 1) and EQ^^. Phylogenetic effect was corrected using 
Felsenstein's independent contrasts'^^. Significance level was set as P< 0.05/21 after 
Bonferroni correction. Species with < 3 ranks were excluded from this analysis to 
minimize statistical errors. EQ for mammals was calculated following the 
allometric formula (Logio(brain weight) — (Logio(body weight) x 0.746-1.253) 
reported by Boddy et al?^, and EQ for non-mammalian amniotes (Logio(brain 
weight) — (Logio(body weight) x 0.55 + 0.0155) were computed following the 
formula reported by Witmer et al?^ 

Analysis of variation in ossification sequence. To examine the rank variation in 
sequence of a particular ossification event, we scaled the rank of each ossification 
event as: 

(r-l)/(r„,„-l) (1) 

in which r is the absolute rank of a given ossification event, and Vax is the total 
number of ranks for each species^^. Therefore, the relative ranks of each species are 
distributed between 0 and 1. This allows removing the differences of maximum 
rank between species resulting from differing levels of sampling resolution between 
species. A similar approach as standardizing the absolute rank r by the maximum 
number of ranks (Vax) has been applied in previous sequence heterochrony 
studies As the ranks are distributed between 1/Vax and 1 with this method, 
the relative ranks of the earliest bone to ossify can vary, depending on Vax- 
However, the method used here circumvents this problem because the relative 
ranks of the earliest event is always be scaled to zero. Nevertheless, some noise 
remains because species with higher have a lower influence on the variance. 
The range in rank variation across species was assessed to examine the variability of 
a particular element in the ossification sequence. As the resolution of the sequence 
can bias the results in this approach, species only with > 3 ranks were included. 

Modularity analysis. Neighbour- joining cluster analysis^^ based on chord 
distance was conducted to identify integration of bone ossification timing. Nodes 
were tested using bootstrapping with 10,000 permutations. Analyses were 
conducted with PAST'^^. Here again, well-resolved species only with > 3 ranks 
were included. Then, three hypothetical module divisions, such as developmental 
modules'^^'^^ (neural-crest-cell bones versus mesoderm bones), ossification mode 
modules (dermal bones versus endochodral bones) and phenotypic variational 
modules'*^'^^, were tested if these could be recovered in neighbour-joining cluster 
analysis. 
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